library(tidyverse)
library(readxl)
library(janitor)
library(broom)
library(lfe)

rm(list= ls())

## 

unemp <- readxl::read_xlsx('data/aa_dat.xlsx') %>%
  mutate(AA.district.id = as.character(gsub("([0-9]+).*$", "\\1", AA_district)))

## Load BT Data 

bt <- readRDS('data/federal_elections_clean.RDS') %>%
  filter(state == 'Bayern') %>%
  mutate(year = as.character(date.id)) %>%
  filter(year %in% c('2013', '2017')) %>%
  left_join(unemp)

## Estimate 2 period DiD for varying RD cutoffs

cutoffs <- seq(.5, 4, .2)

rdd <- lapply(cutoffs, function(bw){
  
  res_right <- bt %>%
    filter(unemp_2015 <= 3.6 + bw & unemp_2015 >= 3.6 - bw) %>%
    felm(right_total_party ~ treated_post + post + treated | 0 |0 | AA.district.id, data = .) %>%
    tidy(conf.int = T) %>%
    filter(term == 'treated_post') %>%
    mutate(outcome = 'Total right vote')
  
  res_left <- bt %>%
    filter(unemp_2015 <= 3.6 + bw & unemp_2015 >= 3.6 - bw) %>%
    felm(left_total_party ~ treated_post + post + treated | 0 |0 | AA.district.id, data = .) %>%
    tidy(conf.int = T) %>%
    filter(term == 'treated_post') %>%
    mutate(outcome = 'Total left vote')
  
  
  return(rbind(res_right, res_left))
  
}) %>% reduce(rbind) %>%
  mutate(bw = rep(cutoffs, each = 2)) 

rdd

## Plot This 

ggplot(rdd, aes(x = bw, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_errorbar(width = 0) + 
  geom_point(shape = 21, fill = 'white') +
  facet_wrap(~ outcome, nrow = 1) +
  theme_bw() + 
  geom_hline(yintercept = 0, linetype = 'dashed') + 
  labs(x = 'Unemployment rate bandwidth',
       y = 'Treatment effect estimate') + 
  scale_x_continuous(breaks = seq(0, 4, 1), 
                     labels = paste(seq(0, 4, 1), '%'))

